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Abstract 

We  have  previously  described  algorithms  for  a  battery  management  system  (BMS)  that  uses  Kalman  filtering  (KF)  techniques  to  estimate  such 
quantities  as:  cell  self-discharge  rate,  state-of-charge,  nominal  capacity,  resistance,  and  others.  Since  the  dynamics  of  electrochemical  cells  are  not 
linear,  we  used  a  nonlinear  extension  to  the  original  KF  called  the  extended  Kalman  filter  (EKF). 

Now,  we  introduce  an  alternative  nonlinear  Kalman  filtering  technique  known  as  “sigma-point  Kalman  filtering”  (SPKF),  which  has  some 
theoretical  advantages  that  manifest  themselves  in  more  accurate  predictions.  The  computational  complexity  of  SPKF  is  of  the  same  order  as  EKF. 
so  the  gains  are  made  at  little  or  no  additional  cost. 

This  paper  is  the  second  in  a  two-part  series.  The  first  paper  explored  the  theoretical  background  to  the  Kalman  filter,  the  extended  Kalman 
filter,  and  the  sigma-point  Kalman  filter.  It  explained  why  the  SPKF  is  often  superior  to  the  EKF  and  applied  SPKF  to  estimate  the  state  of  a 
third-generation  prototype  lithium-ion  polymer  battery  (LiPB)  cell  in  dynamic  conditions,  including  the  state-of-charge  of  the  cell. 

In  this  paper,  we  first  investigate  the  use  of  the  SPKF  method  to  estimate  battery  parameters.  A  numerically  efficient  “square-root  sigma-point 
Kalman  filter”  (SR-SPKF)  is  introduced  for  this  purpose.  Additionally,  we  discuss  two  SPKF-based  methods  for  simultaneous  estimation  of  both 
the  quickly  time- varying  state  and  slowly  time- varying  parameters.  Results  are  presented  for  a  battery  pack  based  on  a  fourth-generation  prototype 
LiPB  cell,  and  some  limitations  of  the  current  approach,  based  on  the  probability  density  functions  of  estimation  error,  are  also  discussed. 

©  2006  Elsevier  B.V.  All  rights  reserved. 

Keywords:  Battery  management  system;  Hybrid  electric  vehicle;  Sigma-point  Kalman  filter;  State-of-charge;  State-of-health;  Lithium-ion  polymer 


1.  Introduction 

This  paper  applies  results  from  the  field  of  study  known  vari¬ 
ously  as  sequential  probabilistic  inference  or  optimal  estimation 
theory  to  advanced  algorithms  for  a  battery  management  system 
(BMS).  This  BMS  is  able  to  estimate  battery  state-of-charge 
(SOC),  instantaneous  available  power,  and  parameters  indica¬ 
tive  of  the  battery  state-of-health  (SOH)  such  as  power  fade  and 
capacity  fade,  and  is  able  to  adapt  to  changing  cell  characteristics 
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over  time  as  the  cells  in  the  battery  pack  age.  The  algorithms  have 
been  successfully  implemented  on  a  lithium-ion  polymer  battery 
(LiPB)  pack  for  hybrid-electric-vehicle  (HEV)  application,2  and 
we  also  expect  them  to  work  well  for  other  battery  chemistries 
and  less  demanding  applications. 

We  have  previously  reported  work  using  extended  Kalman 
filters  (EKF)  to  solve  the  HEV  BMS  algorithm  requirements 
[  1-6].  We  have  since  explored  a  different  form  of  Kalman  filter- 


2  HEV  is  a  particularly  good  benchmark  for  these  algorithms  since  it  imposes 
demanding  requirements  on  a  pack  of  limited  capacity,  resulting  in  cell  elec¬ 
trochemistries  that  are  most  often  far  from  equilibrium.  Therefore,  advanced 
methods  must  be  used  to  estimate  SOC,  SOH,  and  instantaneous  power  in  order 
to  safely,  efficiently  and  aggressively  exploit  the  pack  capabilities. 
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ing  called  sigma-point  Kalman  filters  (SPKF),  and  have  found 
them  to  have  several  important  advantages.  We  introduce  SPKF 
here  in  a  two-part  series,  of  which  this  is  the  second  part.  The 
companion  to  this  paper  [7]  introduces  the  SPKF  and  applies 
it  to  estimating  the  state  of  an  LiPB-based  HEV  battery  cell, 
particularly  its  SOC.  In  this  paper,  we  build  on  the  introduc¬ 
tion  to  first  investigate  the  use  of  the  SPKF  method  to  estimate 
battery  parameters.  A  numerically  efficient  “square-root  sigma- 
point  Kalman  filter”  (SR-SPKF)  is  described  for  this  purpose. 
Additionally,  we  discuss  two  SPKF-based  methods  for  simul¬ 
taneous  estimation  of  both  the  quickly  time-varying  state  and 
slowly  time-varying  parameters.  Applications  to  the  HEV  BMS 
algorithm  requirements  are  outlined,  example  results  given  using 
fourth-generation  prototype  LiPB  cells,  and  conclusions  made. 

We  note  before  continuing  that  there  are  conflicting  objectives 
in  making  this  paper  both  self-contained  and  simultaneously 
reducing  redundant  material  with  respect  to  the  first  paper  in 
this  series.  We  have  perhaps  erred  on  the  side  of  conciseness 
in  several  areas.  For  a  more  in-depth  discussion  of  sequential 
probabilistic  inference  and  the  base-line  SPKF  algorithm,  the 
reader  is  referred  to  Ref.  [7], 

2.  Summary  of  sequential  probabilistic  inference 

Very  generally,  any  causal  dynamic  system  (e.g.,  a  battery 
cell)  generates  its  outputs  as  some  function  of  its  past  and  present 
inputs.  Often,  we  can  define  a  state  vector  for  the  system  whose 
values  together  summarize  the  effect  of  all  past  inputs.  Present 
system  output  is  a  function  of  present  input  and  present  state 
only;  past  input  values  need  not  be  stored.  The  system’s  pa¬ 
rameter  vector  comprises  all  quasi-static  numeric  quantities  that 
describe  how  the  system  state  evolves  and  how  the  system  out¬ 
put  may  be  computed.  The  state-vector  quantities  change  on  a 
relatively  rapid  time  scale,  and  the  parameter-vector  quantities 
change  on  a  relatively  long  time  scale  (or,  not  at  all). 

We  assume  that  the  electrochemical  cell  under  consideration 
may  be  modeled  using  a  discrete-time  state-space  model  of  the 
form 

Xk  =  fix, t_i,  Uk—\ ,  tVk—l ,  k  -  1)  (1) 

yk  =  h(xk,  Uk,Vk,  k).  (2) 

Here,  xk  e  R"  is  the  system  state  vector  at  time  index  k.  and  Eq. 
(1 )  is  called  the  “state  equation”  or  “process  equation”  and  cap¬ 
tures  the  evolving  system  dynamics.  The  known/deterministic 
input  to  the  system  is  uk  e  Rp,  and  wk  e  M"  is  stochastic  “pro¬ 
cess  noise”  or  “disturbance”  that  models  some  unmeasured  input 
which  affects  the  state  of  the  system.  The  output  of  the  sys¬ 
tem  is  yk  e  M'" ,  computed  by  the  “output  equation”  (2)  as  a 
function  of  the  states,  input,  and  Vk  e  Rm,  which  models  “sen¬ 
sor  noise”  that  affects  the  measurement  of  the  system  output 
in  a  memoryless  way,  but  does  not  affect  the  system  state. 
f(xk- 1,  Uk- 1,  Wk- 1,  k  —  1)  is  a  (possibly  nonlinear)  state  tran¬ 
sition  function  and  g(xk,  Uk ,  Vk,  k)  is  a  (possibly  nonlinear)  mea¬ 
surement  function. 

In  the  companion  to  this  paper  we  introduced  the  relevant  the¬ 
ory  describing  the  optimal  solution  to  state  estimation  -  based  on 


Table  1 

Summary  of  the  general  sequential  probabilistic  inference  solution 

General  state-space  model 

Xk  =  fi*k- 1,  uk- 1,  wk-i,k  -  1) 
yk  =  H*k,  uk,  Vk,  k), 

where  wk  and  v k  are  independent,  Gaussian  noise  processes  of  covariance 
matrices  Ew  and  Sv,  respectively 

Definitions:  let 

xf=xk-x~,  %  =  yk-  % 

Initialization:  for  k  =  0,  set 
x o  =  EHol 

r+0  =  E[Uo-x+)(A-0-i-+)T] 

Computation:  for  k  =  1,2,...  compute 

State  estimate  time  update:  xf  =  E[/(x^_i ,  uk-\ ,  wk-i,  k  —  l)|Yjt_i] 
Error  covariance  time  update:  £fk  =  E[(^)(jc^)T] 

Output  estimate:  yk  =  E [h(xk,  uk,  vk ,  k )  |  Y^_i] 

Estimator  gain  matrix:  Lk  =  E[(ip(Ji)T](E[(yH(»)T])-1 
State  estimate  measurement  update:  xf  =  +  Lk(yk  —  yk) 

Error  covariance  measurement  update:  -LV =  27“  —  LkSy<kLj 


a  body  of  knowledge  known  as  “sequential  probabilistic  infer¬ 
ence”  -  and  showed  how  various  approximations  and  assump¬ 
tions  resulted  in  the  Kalman  filter,  the  extended  Kalman  filter, 
and  the  sigma-point  Kalman  filter  [7],  Here,  we  summarize  the 
key  points  required  to  understand  parameter  estimation  and  si¬ 
multaneous  state  and  parameter  estimation.  For  greater  detail, 
the  reader  should  consult  Ref.  [7], 

The  general  solution  from  which  all  the  KF  variants  are  de¬ 
rived  comprises  two  major  stages  per  measurement  interval. 
First,  it  predicts  the  present  state  (the  “time  update”)  given  all 
prior  information;  second,  it  corrects  the  prediction  using  the 
current  measurement  (the  “measurement  update”).  The  KFs  not 
only  estimate  the  state,  but  also  the  state  error  covariance  matrix 
to  provide  an  ongoing  uncertainty  estimate  on  the  state  estima¬ 
tion  error. 

Table  1  summarizes  the  general  solution.  In  the  notation  we 
use,  the  decoration  “circumflex”  indicates  an  estimated  quantity 
(e.g.,  x  indicates  an  estimate  of  the  true  quantity  x ).  A  super¬ 
script  ”  indicates  an  a  priori  estimate  (i.e.,  a  prediction  of 
a  quantity’s  present  value  based  on  past  data)  and  a  superscript 
“+”  indicates  an  a  posteriori  estimate  (e.g.,  xj  is  the  estimate  of 
true  quantity  x  at  time  index  k  based  on  all  measurements  taken 
up  to  and  including  time  k).  The  decoration  “tilde”  indicates  the 
error  of  an  estimated  quantity. 

Y*  =  [yk,  yk- 1,  ■  •  • ,  vo}- 

The  symbol  Exy  =  E[ryT]  indicates  the  auto-  or  cross¬ 
correlation  of  the  variables  in  its  subscript.  (Note  that  often  these 
variables  are  zero-mean,  so  the  correlations  are  identical  to  co- 
variances.)  Also,  for  brevity  of  notation,  we  often  use  Xx  to 
indicate  the  same  quantity  as  Uxx. 

The  algorithm  first  initializes  the  filter  and  iteratively  ex¬ 
ecutes  six  steps  each  measurement  interval.  First,  the  present 
state  value  is  predicted  using  prior  data.  Second,  the  covariance 
of  this  state  estimate  error  is  updated.  Third,  the  present  cell  out¬ 
put  is  predicted.  Fourth,  the  Kalman  gain  matrix  is  computed. 
Fifth,  the  actual  cell  output  is  measured  and  compared  to  the 
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Table  2 

Weighting  constants  for  two  sigma-point  methods 


Y 


(m) 

*0 


(m) 

4 


(c) 

v 


a 


(c) 

i 


UKF 

CDKF  h 


k  1 

L  +  k  2(L  +  k) 
h2-  L  1 
h- 2  h2 


y-rr  +  a  -  “2  +  F> 

L  t  A 

h2  -  L 
~h2~ 


1 

2 (L  +  k) 

1 

lii2 


k  =  a2(L  +  k)  —  L  is  a  scaling  parameter,  with  (1 0  2  <  a  <  1).  Note  that  this 
a  is  different  from  a'"'1  and  <yici .  at  is  either  0  or  3  —  L .  /3  incorporates  prior  infer- 
mation.  For  Gaussian  RVs,  /3  =  2.  h  may  take  any  positive  value.  For  Gaussian 
RVs,  h  =  V3. 


estimate.  The  output  estimate  error  is  weighted  by  the  Kalman 
gain  and  used  to  update  the  state  estimate.  Finally,  the  state  esti¬ 
mate  error  covariance  is  updated.  The  output  of  the  filter  at  each 
time  instant  is  the  state  estimate  from  step  5  and  the  error  co- 
variance  from  step  6  (used  to  compute  error  bounds  on  the  state 
estimate). 

3.  Summary  of  sigma-point  Kalman  filters 

For  nonlinear  systems,  a  closed-form  solution  or  even  an  al¬ 
gorithm  to  exactly  implement  the  general  probabilistic  inference 
solution  in  Table  1  is  not  available.  The  EKF  is  one  approach 
to  approximating  the  solution  using  a  first-order  linearization  of 
the  system  dynamics,  which  is  based  on  some  questionable  as¬ 
sumptions.  Sigma-point  Kalman  filters  are  an  alternate  approach 
to  generalizing  the  Kalman  filter  to  state  estimation  for  nonlin¬ 
ear  systems.  SPKFs  rely  on  numeric  approximations  rather  than 
analytic  approximations  of  the  EKF.  Each  time  step,  a  set  of 
points  ( sigma  points )  is  chosen  so  that  the  (possibly  weighted) 


mean  and  covariance  of  the  points  exactly  matches  the  mean  and 
covariance  of  the  a  priori  random  variable.  These  points  are  then 
passed  through  the  nonlinear  function,  resulting  in  a  transformed 
cloud  of  points.  The  a  posteriori  mean  and  covariance  that  are 
sought  are  then  approximated  by  the  mean  and  covariance  of  this 
cloud.  Note  that  the  sigma  points  comprise  a  fixed  small  num¬ 
ber  of  vectors  that  are  calculated  deterministically — not  like  the 
Monte  Carlo  or  particle  filter  methods.  This  method  is  used  to 
approximate  state  and  output  estimates,  and  their  covariances. 

Specifically,  if  the  input  random  vector  x  has  dimension  L, 
mean  x,  and  covariance  X*,  then  p  +  1  =  2L  +  1  sigma  points 
are  generated  as  the  set 

X  =  [x,  x  +  Ysf^x,  x  -  yfz*}, 

with  elements  of  X  indexed  from  0  to  p.  and  where  the  matrix 
square-root  R  —  computes  a  result  such  that  X  =  RR[ . 
Usually,  the  efficient  Cholesky  decomposition  [8,9]  is  used,  re¬ 
sulting  in  lower-triangular  R.  The  reader  can  verify  that  the 
weighted  mean  and  covariance  of  X  equal  the  original  mean  and 
covariance  of  random  vector  x  for  a  specific  set  of  { y,  o'""*,  alc>) 
if  we  define  the  weighted  mean  as  x  —  Ef=o  c^fXi,  the 
weighted  covariance  as  Xj  =  Ef=o  af(Xi  —  x)(Xj  —  x)T,  Xj 
as  the  zth  element  of  X,  and  both  a'""  and  af  as  real  scalars  with 
the  necessary  (but  not  sufficient)  conditions  that  E?=o  0!;m,  =  1 
and  ZLo  —  1  ■  The  various  sigma-point  methods  differ  only 
in  the  choices  taken  for  these  weighting  constants.  Values  for  the 
two  most  common  methods  -  the  unscented  Kalman  filter  (UKF) 
[10-15]  and  the  central  difference  Kalman  filter  (CDKF)  [16- 
18]  -  are  summarized  in  Table  2.  In  this  work,  we  use  the  CDKF 
parameters. 


Table  3 

Summary  of  the  nonlinear  sigma-point  Kalman  filter 


Nonlinear  state-space  model 

*k  =  f(Xk- 1,  «i-l,  U>k-u  k  -  1) 
yk  =  h(xk,  Mi,  t>t,  k) 

where  w k  and  Vk  are  independent,  Gaussian  noise  processes  of  covariance  matrices  Uw  and  Uv,  respectively 


Definitions:  let 

4  =  [■**.  wl  ^  x“k  =  K**)T.  (*i)T.  (Ti)TlT, 

Initialization:  for  k  =  0,  set 
xf  =  E[*ol 

S+0  =  E[(*o-*J)Uo--tf)T] 

Computation:  for  k  =  1,2,...  compute 
State  estimate  time  update 


p  =  2  x  dim(jtjp 


Error  covariance  time  update 
Output  estimate 


Estimator  gain  matrix 


State  estimate  measurement  update 
Error  covariance  measurement  update 


xf+  =  E[jrg]  =  [(*0  )T,  w,  u]1 

Kt  =  EK*o  - *o +X*o  “ *o+)Tl  =  diag(r+0, Sw,  Sv) 


xk- 1 

=  {*t+i>*i-i  + 

xk- 1  -  Y 

xx,~ 

k,i 

=  f<*l 

’Xk-U’k- 

-  1) 

K  = 

V''  a(m)Xx’~ 
Z^/=0  01  i  *  k,i 

Zf,k 

H 

* 

.3 

W 

ii 

-*i)T 

yk,i  = 

=  h(Xxk-~ ,  uk,  Xf_ 

\,vk) 

%  = 

ELo  4m)y« 

Zy.k 

=  Ef=o“ic)(Xi.,' 

-  5 4)04,;  - 

Vi)7 

E7.  , 

xy,k 

r-vxTi, 

i  -  yk)T 

Lk  = 

Sxy,kSy,k 

=  *k  +  Lk(yk  -  h) 
stk  =  Sx,k  -  LkEy,kLj 
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To  use  SPKF  in  an  estimation  problem,  we  first  define  an 
augmented  random  vector  xa  that  combines  the  randomness  of 
the  state,  process  noise,  and  sensor  noise.  This  augmented  vector 
is  used  in  the  estimation  process.  The  math  behind  the  SPKF  is 
described  in  Ref.  [7]  and  the  overall  solution  is  summarized  in 
Table  3. 


4.  Computationally  efficient  square-root  sigma-point 
Kalman  filters 


Sigma-point  Kalman  filters  may  be  used  directly  for  state 
estimation  and  we  have  shown  that  they  produce  better  state 
estimates  and  much  better  covariance  estimates  than  EKF  [7]. 
The  computational  complexity  is  (7(Z.3),  which  is  of  equivalent 
complexity  to  EKF  state  estimation,  where  L  is  the  dimension 
of  the  augmented  state.  They  may  also  be  used  for  parameter 
estimation,  as  will  be  described  in  the  sequel,  but  the  compu¬ 
tational  complexity  remains  0(L 3),  whereas  the  corresponding 
EKF  method  has  complexity  0(L2). 

The  bottleneck  in  the  SPKF  algorithm  is  the  computation  of 
the  matrix  square-root  SkSj  =  Ek  each  time  step,  which  has 
computational  complexity  0(Li / 6)  using  a  Cholesky  factor¬ 
ization.  A  variant  of  the  SPKF,  to  which  we  will  refer  as  the 
square-root  SPKF,  propagates  Sk  directly  without  needing  to 
re-factor  each  time  step  [15,19,20].  This  approach  has  several 
advantages:  there  are  improved  numeric  properties  as  the  covari¬ 
ances  are  guaranteed  to  be  positive  semi-definite,  and  although 
the  state  estimation  update  is  still  0(L 3),  the  parameter  update 
may  now  be  done  in  0(L2).  Therefore,  for  the  same  computa¬ 
tional  complexity  of  EKF,  better  results  are  obtained. 

Three  important  linear-algebra  techniques  are  required  to  im¬ 
plement  SR-SPKF  [8,9]: 


•  QR  decomposition :  The  QR  decomposition  algorithm  com¬ 
putes  two  factors  Q  e  RNxN  and  R  e  RLxN  for  a  matrix 
A  eS.LxN  such  that  A  =  QR,  Q  is  orthogonal,  R  is  upper- 
triangular,  and  N  >  L.  The  property  of  the  QR  factorization 
that  is  important  here  is  that  R  is  related  to  the  Cholesky  fac¬ 
tor  we  wish  to  find.  Specifically,  if  R  e  is  the  upper- 

triangular  portion  of  R,  then  A'1  is  the  Cholesky  factor  of 
E  =  ATA.  That  is,  if  R  —  qr(AT)T,  where  qr(-)  performs  the 
QR  decomposition  and  returns  the  upper-triangular  portion 
of  R  only,  then  R  —  chol(AAT).3 

This  information  may  be  used  by  noticing  that  step  2  of 
the  SPKF  method  computes 


ru  =  E“/cW’ 
1=0 


\X 


k,i 


hm:7 


rx  1T 
xk,i>  > 


which  may  also  be  written  as  E-  k  =  AAT,  where 

/ (c)  X  — 

A  =  [V<'(^\0:p)-%“  )].  A  similar  computation  is 
performed  in  step  4.  Rather  than  computing  AAT  and  then 


3  Some  implementational  care  is  advised  as  the  QR  decomposition  and 
Cholesky  factorizations  are  not  defined  identically  by  all  numeric  packages. 
For  example,  chol()  or  qr()  may  return  the  transpose  of  the  desired  result. 


later  computing  the  Cholesky  factor  thereof,  we  can  instead 
compute  the  QR  decomposition  of  AT,  and  propagate  the  R 
component.  The  computational  complexity  of  the  QR  decom¬ 
position  is  0(NL2),  whereas  the  complexity  of  the  Cholesky 
factor  is  0{L2 / 6)  plus  0(NL2)  to  first  compute  AAT. 

•  Cholesky  downdating :  Note  that  the  previous  method 
fails  whenever  a)  is  negative,  because  the  square-root 
involved  will  not  have  a  real  value.  In  particular,  aX  is 
frequently  negative,  and  the  final  step  of  SPKF  computes 
E- !",.  =  Ekk  —  LkEykLj.  Using  the  latter  case  as  an 
example,  we  cannot  simply  append  a  column  to  A  with 
value  yj— Ey  kLk  because  of  the  negative  sign  inside  the 
square-root.  The  solution  to  this  problem  is  the  Cholesky 
downdating  procedure,  which  takes  as  input  the  Cholesky 
factor  to  A  At  and  Q  Ey^Lk,  and  computes  the  updated  factor 
for  AAt  —  LkEy^Lj.  If  Eyik  is  not  a  scalar,  downdating  is 
done  using  each  of  its  columns.  This  algorithm  is  only  0{L2). 

•  Backsubstitution :  Finally,  the  solution  to  the  estimator  gain 
matrix  computation  is  LkEyk  =  EE  k,  which  may  be  written 

as  Lk  =  Z7y,k(S~y,ksJ.,k)~1’  or  alternately  as  Lk(,SyikSjk )  = 
EE  k.  If  yk  is  not  a  scalar,  this  equation  may  often  be 
computed  most  efficiently  via  back-substitution  in  two  steps. 
First,  (z)Sj,k  =  EEyk  is  found,  and  then  LkS-y^  =  Z  is 
solved.  Since  Sy, k  is  already  triangular,  no  matrix  inversion 
need  be  done.  Backsubstitution  has  complexity  0{N2 / 2). 


To  use  SR-SPKF  in  an  estimation  problem,  we  again  define  an 
augmented  random  vector  xa  that  combines  the  randomness  of 
the  state,  process  noise,  and  sensor  noise.  This  augmented  vector 
is  used  in  the  estimation  process  as  described  below. 

SR-SPKF  step  1:  state  estimate  time  update.  As  with 
SPKF,  each  measurement  interval,  the  state  estimate  time 
update  is  computed  by  first  forming  the  augmented  a  pos¬ 
teriori  state  estimate  vector  for  the  previous  time  interval: 
xakf  —  [(AjLi)T,  w,  i)]T.  With  SR-SPKF,  the  square-root  aug¬ 


mented  a  posteriori  covariance  estimate  is  computed:  S?k_  \  — 
diag(S^T_j,  Sw,  Sv).  These  factors  are  used  to  generate  the 
P+  1  sigma  points:  X°kf  =  {xkf,xakf  +  yS~f_ vxakf  - 
ySk'fi  }  •  From  the  augmented  sigma  points,  the  p  +  1  vectors 


comprising  the  state  portion  Xk  \  and  the  p  +  1  vectors  com 


prising  the  process  noise  portion  Xkf  are  extracted.  The  pro¬ 
cess  equation  is  evaluated  using  all  pairs  of  Xkf  .  and  Xkf  . 
(where  the  subscript  i  denotes  that  the  zth  element  is  being  ex¬ 
tracted  from  the  original  set),  yielding  the  a  priori  sigma  points 
Xk'~  for  time  step  k.  Finally,  the  a  priori  state  estimate  is  com¬ 
puted  as  xk  =  £f=o  u^Xl’j. 

SR-SPKF  step  2:  error  covariance  time  update.  Using  the  a 
priori  sigma  points  from  step  1,  the  square -root  a  priori  covari¬ 
ance  estimate  is  computed  as 


SR-SPKF  step  3:  estimate  system  output  yk .  The  system  out¬ 
put  is  estimated  by  evaluating  the  model  output  equation  using 
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Table  4 

Summary  of  the  square-root  sigma-point  Kalman  filter  for  state  estimation 

Nonlinear  state-space  model 
rk  =  f(*k- 1,  Uk-i,  wk-\,k  -  1) 
yk  =  h(xk,  uk,  vk,  k) 

where  Wk  and  Vk  are  independent,  zero-mean,  Gaussian 

noise  processes  of  covariance  matrices  Sw  and  Zv,  respectively 

Definitions:  let 

4  =  [*?.  toj.  vjf,  X%  =  1  (Xl)\  (*»)T,  (Xm 

T,  p  =  2  x  dim(jp 

Initialization:  for  k  =  0,  set 

X, J  =  E[*ol 

sto  =  ch°l{EK*o  -  *o  X-to  -  i0+)Tl) 

*o'+  =  EP'ol  =  K*o  E  °-  °1T 

=  chol{E[(x“  -  x“-+)0"  -  *“'+)T]}  =  diag(Sf( 

Computation:  for  k  =  1,2,...  compute 

State  estimate  time  update 

K-1  =  +  y%tv%- L  -  yS&- 1> 

xi,7  =  “*-!■  xk-u'  k~X) 

a*-  =  £  L-Lk; 

Error  covariance  time  update 

^*  =  qr{[ 

Output  estimate 

yk,i  =  h(Xxk’~,  Uk,  xvk’+  .,  k) 

s*=£?=o  *T°yu 

Estimator  gain  matrix 

Sy.k  =  qr{[\/apQ4,(0:p)  -  »)T11T 

=  EL  -  w  -  »)T 

Lk  =  E~~  k(Sy  k^yj)  1  (s°lve(i  by  backsubstitution) 

State  estimate  measurement  update 

Error  covariance  measurement  update 

K  =  K  +  -  h) 

slk  =  downdate), SJj.,  LkSy,k 1 
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the  sigma  points  describing  the  spread  in  the  state  and  noise  vec¬ 
tors. First,  wecompute  the  points  34,/  =  h(Xk’P,  uk,  X'k- 1  i ’  ^)- 

The  output  estimate  is  then  %  —  El’=o  34,/- 

SR-SPKF  step  4:  estimator  gain  matrix  Lj..  To  compute  the 
estimator  gain  matrix,  we  must  first  compute  the  required  co- 
variance  and  square-root  covariance  matrices. 
r  r  , -  -i  >  T 


Sy ,k  =  qr  | 

Zry,k  =  E 


p 

1=0 


otf\yk,(0:P)  -  h)T  j 
ctf\xxk’7  -  xk){yu  -  %). 


Then,  we  simply  compute  Lk  =  k^yl”  s°lve(i  by  backsub- 
stitution. 

SR-SPKF  step  6:  error  covariance  measurement  update.  The 
final  step  computes  the  square-root  form  of  Sf  —  Skk  — 

LkUy.kLj  as  Sfk  —  downdate  j Sk  k,  L*S-y,£  j.  The  SR-SPKF 
solution  for  state  estimation  is  summarized  in  Table  4,  and  re¬ 
sults  for  both  SPKF  and  SR-SPKF  for  SOC  estimation  are  given 
in  Section  8.3. 


5.  Parameter  estimation  using  SR-SPKF 

The  various  methods  for  state  estimation  presented  so  far  have 
assumed  a  known  system  model  in  the  form  of  Eqs.  (1)  and  (2). 
These  equations  will  generally  involve  numeric  values  in  their 
computations.  Some  of  these  values  might  be  intrinsic  constants, 
but  others  might  be  factors  determined  by  the  electrochemistry 
or  construction  of  a  particular  cell.  We  refer  to  these  latter  factors 
as  the  “parameters”  of  the  cell  model. 

To  use  the  enhanced  self-correcting  (ESC)  cell  model  as  an 
example  (cf.  Section  7),  the  parameters  comprise  the  following: 
the  Coulombic  efficiency  rj,  the  total  capacity  C,  the  filter  poles 


a i ,  . . . ,  allf,  the  filter  weighting  factors  gi, . . . ,  gnf_i ,  the  cell 
discharge  and  charge  resistances  R+  and  R~ ,  the  hysteresis  rate 
constant  y,  and  the  maximum  level  of  hysteresis  M.  Combined, 
they  are 

0  =  [p,  C,a i, . . .  ,an/,  gi, . . . ,  gn f-l ,  R+ ,  R~ ,  y,  M]J . 

We  have  previously  shown  how  to  estimate  a  system’s  state 
given  a  known  model  and  noisy  measurements.  We  now  show 
how  to  estimate  a  system’s  parameters  given  a  known  state  and 
clean  measurements.  We  assume  that  there  is  a  true  value  for  9 
that  describes  the  cell  under  consideration,  and  wish  to  adapt  an 
estimate  6  to  converge  to  the  true  value.  We  begin  by  proposing  a 
state-space  model  for  the  “dynamics”  of  the  system  parameters. 


Ok  =  Ok- 1  +  rk- 1 

dk  =  h(xk,  uk,  Vk,  Ok,  k )  +  ek- 

The  first  equation  states  that  the  parameters  are  essentially  con¬ 
stant,  but  that  they  may  change  slowly  over  time  by  some  driving 
process,  modeled  by  a  process  rk  of  small  fictitious  “noise”.  The 
output  equation  for  the  state-space  model  of  true  parameter  dy¬ 
namics  is  the  estimate  of  cell  output  based  on  the  previous  state, 
the  measured  inputs,  and  parameter  estimates,  added  to  some 
estimation  error  e*. 

With  this  system  defined,  we  can  apply  SR-SPKF  to  estimate 
the  parameters.  The  SR-SPKF  solution  for  parameter  estimation 
is  summarized  in  Table  5.  The  dominant  differences  between 
state  and  parameter  estimation  are:  in  parameter  estimation  the 
noises  are  assumed  to  contribute  linearly  to  the  state  and  mea¬ 
surement  equations,  and  an  approximate  process  noise  update  is 
performed  (see  [  15]  for  details). 
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Table  5 

Summary  of  the  square-root  sigma-point  Kalman  filter  for  parameter  estimation 


Nonlinear  state-space  model 
Ok  =  Ok- 1  +  fk-l 
dk  =  h(xk,  uk,  vk,  0k,  k )  +  ek 

where  rk  and  ek  are  independent,  zero-mean,  Gaussian  noise  processes  of  covariance  matrices  Ur  and  Se,  respectively 


Definition:  let 

Dr,k- 1  =  -diagjS^  j)  +  f  diagjS^  j}2  +  diag{ 1 }  p  =  2x  dim(ft) 

Initialization:  for  k  =  0,  set 
e+  =  E[0O]  s+0  =  Choi  {  E[(flo  -  Sjjm  -  e+ f] } 

Computation:  for  k  =  1,2,...  compute 
Parameter  estimate  time  update 
Error  covariance  time  update 

Output  estimate 


Estimator  gain  matrix 


Parameter  estimate  measurement  update 
Error  covariance  measurement  update 


«*“  =  el 
sz,  =  s ; 


e,k- 1 


+  A 


wk  =  {e-,e^  +  YS7re--Ysik} 

Vk,i  =  h(xk,  uk,vk,  Wi,i,  k) 

4  =  EfAmA 


Sd  k  =  qr{[ Af  (2?*,(0:p)  -  dk)  v^]T)T 

sh,k  =  ELo  “f A,-  -  K)(Vkj  -  dk? 

Lk  =  S—  (S?  Sj  k)-1  (solved  by  backsubstitution) 
=  0k  +  Lk(dk  —  dk) 

Stk  =  downdatefEr  LkS~d.k) 


6.  Joint  and  dual  sigma-point  filtering 

We  have  now  shown  how  to  estimate  the  state  of  a  system 
given  a  known  model  and  noisy  measurements  and  how  to  esti¬ 
mate  the  parameters  of  the  system  given  a  known  state  and  clean 
measurements.  We  proceed  to  give  two  methods  whereby  one 
can  simultaneously  estimate  both  the  state  and  parameters  of  a 
system  given  noisy  measurements. 

The  various  methods  of  state  estimation  presented  so  far  have 
assumed  a  constant  system  model.  However,  when  applying 
these  procedures  to  estimate  battery  SOC,  for  example,  we  en¬ 
counter  a  possible  source  of  error:  not  all  cells  are  created  equal. 
Most  of  the  research  reported  here  has  been  conducted  using 
prototype  high-power  LiPB  cells,  constructed  by  hand.  In  these 
cells,  there  is  a  great  deal  of  variability  in  resistance,  capacity, 
time  constants,  and  so  forth.  Even  when  mass-produced,  how¬ 
ever,  there  is  some  cell-to-cell  variation,  which  only  increases 
as  the  cells  age,  both  in  accumulated  cycles  and  in  calendar  life. 

Some  of  the  critical  parameters,  such  as  cell  resistance  and 
capacity,  directly  limit  the  pack  performance  through  “power 
fade”  and  “capacity  fade”,  so  are  indicators  of  the  battery  state- 
of-health.  It  is  important  to  be  able  to  estimate  these  and  other 
parameters  to:  (1)  maintain  an  accurate  model  for  SOC  estima¬ 
tion,  (2)  understand  the  present  battery  state-of-health,  and  (3) 
predict  remaining  service  life. 

Keeping  in  mind  the  previous  discussion  on  estimating  SOC, 
it  is  apparent  that  the  quantities  descriptive  of  the  present  battery 
pack  condition  exist  on  two  time  scales.  Some  change  rapidly, 
such  as  SOC,  which  can  traverse  its  entire  range  within  minutes. 
Others  may  change  very  slowly,  such  as  cell  capacity,  which 
might  change  as  little  as  20%  in  a  decade  or  more  of  regular 
use.  The  quantities  that  tend  to  change  quickly  comprise  the 
state  of  the  system,  and  the  quantities  that  tend  to  change  slowly 
comprise  the  time-varying  parameters  of  the  system. 


Any  of  the  KF  methods  used  to  estimate  SOC  might  be 
adapted  to  concurrently  estimate  both  the  state  and  the  slowly 
time-varying  cell  parameters  by  combining  the  cell  model  state 
vector  with  the  model  parameters  and  simultaneously  estimat¬ 
ing  the  values  of  this  augmented  state  vector.  This  method  is 
called  joint  estimation.  It  has  the  disadvantages  of  large  matrix 
operations  due  to  the  high  dimensionality  of  the  resulting  aug¬ 
mented  model  and  potentially  poor  numeric  conditioning  due  to 
the  vastly  different  time  scales  of  the  states  (including  parame¬ 
ters)  in  the  augmented  state  vector.  However,  it  is  quite  straight¬ 
forward  to  implement.  We  first  combine  the  state  and  parameter 
vectors  to  form  an  augmented  state  such  that  the  dynamics  may 
be  represented  by 


Xk 

f(xk-i,uk-i,dk-i,  Wk-I,k  -  1) 

A 

dk- 1  +  n~  i 

yk  =  h(xk,  Uk ,  Vk,  k). 


Note  that  to  simplify  notation,  we  will  refer  to  the  vector  com¬ 
prising  both  the  present  state  and  the  present  parameters  as  X*, 
the  vector  comprising  the  present  process  noise  and  present  pa¬ 
rameter  noise  as  W*,  and  the  equation  combining  the  dynamics 
of  the  state  and  the  dynamics  of  the  parameters  as  T.  This  allows 
us  to  write 

Xk  =  f(Xn,  uk-\,  Wjfc_i,  k  -  1) 

yk  =  h{Xj t,  uk,  Vk,  k). 

With  the  augmented  model  of  the  system  state  dynamics  and 
parameter  dynamics  defined,  we  apply  the  SPKF  method.  Table 
6  gives  a  listing  of  the  steps,  which  correspond  directly  to 
the  state  estimation  SPKF  steps,  except  with  larger  matrix 
operations. 

The  second  method,  dual  estimation,  also  uses  a  Kalman  fil¬ 
tering  approach  to  estimate  both  state  and  parameter  values,  but 
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Table  6 

Summary  of  the  nonlinear  sigma-point  Kalman  filter  for  joint  estimation 
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Xk 

f(xk-i,  Uk- 1,  Wk-1, 0k-l,  k  -  1) 

ok 

8k~l  +  rk- 1 

State-space  model 

i  hi.  i  mi.  i  f)i.  i  k  —  1 

X,t  =  W*_i,  k  —  1) 

yt  =  h(Xk,  uk,  vk ,  k) 

yk  =  h(xk,  uk,  vk,  9k,  k) 

where  Wk,rk,  and  Vk  are  independent,  Gaussian  noise  processes  of  covariance  matrices  Zw,  Zr,  and  Zv,  respectively.  For 
brevity,  we  let  X*  =  [xj,  #J]T,  W k  =  [wj ,  rJ]T  and  Zyy  =  diag(X'U),  Er) 

Definitions:  let 

X“  =  [XT  WT  vjf,  Xak  =  [(Af)T,  (A’W)T,  (Ar")T]T,  p  =  2  x  dim(X“) 


Initialization:  for  k  =  0,  set 

X+  =  E[X0] 

X+,  =  E[(Xo  -  X+)(X0  -  X+)T] 
Computation:  for  k  =  1,2,...  compute 
State  estimate  time  update 


Error  covariance  time  update 
Output  estimate 

Estimator  gain  matrix 


State  estimate  measurement  update 
Error  covariance  measurement  update 


x“-+ 

=  E[Xg] 

=  1(X+)T, 

W,  D]t 

^x.o 

=  E[(Xg 

—  Xq +)(Xq  —  Xq +)T]  = 

:dia^i,<r 

ya>+ 

Ak- 1 

=  {* 

ita'+  -U 
■V  Afc-1  + 

“  Y\[' 

<r 

=  nx\ 

-l.i'  uk-l' 

X^ti.k-l) 

x-  = 

-EL  0“ 

(m)  yX,— 

/  Ak,i 

si.k 

=  EL> 

4\xf:r 

-w*r- 

-x-)T 

yk.i  -- 

=  *(*«■ 

1-+ 

.».*> 

%  = 

Zy.k 

= 

Z-/i=o 

af\yt,i  - 

yk)(yk,i  -  %) 

r 

=  TP 

■  -  x-)(w,i  - 

»>T 

Lk  = 

z~  z: 

Xy,k  } 

-1 

\k 

•Ew.  Ev) 


^ k  —  +  Lk(yk  —  yk) 
^k  =  ^x,k  ~  Lk£y,kLj 


uses  separate  Kalman  filters  for  state  estimation  and  parameter 
estimation.  The  computational  complexity  is  smaller  and  the  ma¬ 
trix  operations  may  be  numerically  better  conditioned.  However, 
by  decoupling  state  from  parameters,  any  cross-correlations 
between  changes  are  lost,  leading  to  potentially  poorer 
accuracy. 

The  mathematical  model  of  cell  dynamics  again  explicitly 
includes  the  parameters  as  the  vector  Ok : 

Xk  =  f(xk-l,  Uk-l,  Wk-i,0k-i,k  -  1) 
yk  =  h(xi t,  uk,  Vk,  9k-i ,k). 

Non-time-varying  numeric  values  required  by  the  model  may 
be  embedded  within  /(•)  and  li(-),  and  are  not  included  in  6k. 
We  also  slightly  revise  the  mathematical  model  of  parameter 
dynamics  to  explicitly  include  the  effect  of  the  state  equation. 

Ok  —  0k-i  +  rk- 1 

dk  =  h(f(xk- 1,  Uk-i,  Wk-i,0k-i,k  -  1),  uk,  Vk,  9k-i,k)+ek- 

With  these  two  systems  defined,  we  can  apply  the  standard  pro¬ 
cedure  of  dual  extended  Kalman  filtering  [6,21],  or  generalize 
the  procedure  to  other  forms  of  Kalman  filtering.  Dual  sigma- 
point  Kalman  filtering  is  outlined  in  Table  7,  and  comprises  two 
carefully  integrated  SPKFs.  The  algorithm  is  initialized  with  the 
best  guess  of  the  true  parameters  =  E[0q].  and  with  the  best 
guess  of  the  cell  state  xj"  =  E[xq].  The  estimation  error  covari¬ 
ance  matrices  are  also  initialized.  The  algorithm  may  be  adapted 
to  use  an  SR-SPKF  for  the  state-  and/or  the  parameter-estimation 
filter  in  a  straightforward  manner. 

The  dual  sigma-point  Kalman  filter  can  be  viewed  by  drawing 
a  block  diagram,  as  in  Fig.  1.  We  see  that  the  process  essentially 


comprises  two  sigma-point  Kalman  filters  running  in  parallel  - 
one  adapting  the  state  and  one  adapting  parameters  -  with  some 
information  exchange  between  the  filters. 

6.1.  Convergence 

The  dual/joint  extended  Kalman  filters  will  adapt  %  and  6k  so 
that  the  model  input-output  relationship  matches  the  cell  input- 
output  data  as  closely  as  possible.  There  is  no  built-in  guarantee 
that  the  state  or  parameters  of  the  model  converge  to  anything 
with  physical  meaning.  We  take  special  steps  to  ensure  that  this 
occurs. 

A  very  crude  cell  model  may  be  used,  combined  with  the 
dual/joint  SPKF,  to  ensure  convergence  of  the  SOC  state.  Specif- 


Fig.  1.  Diagram  of  dual  estimation  method.  Solid  lines  represent  state-  and 
parameter-vector  signal  flow,  and  dashed  gray  lines  represent  error  covariance 
matrix  signal  flow. 
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Table  7 

Summary  of  the  sigma-point  Kalman  filter  for  dual  state  and  parameter  estimation 


Nonlinear  state-space  models 

xk  =  K*_1,  Wk-1,  0k-i,  k  -  1)  and  ek  =  0*_i  +  rk- 1 

yk  =  h(xk,  uk,  vk,  9k ,  k)  dk  =  h(f(xk-i,uk-i,  ui*-i,  ftt-i -k  -  1),  uk,  vk,  6k-i,k)  +  ek 

where  Wk,  Vk,  rk  and  ek  are  independent,  Gaussian  noise  processes  of  covariance  matrices  Zw,  Uv,  Er  and  Ue,  respectively 


Definitions 

Xak  =  [. xj ,  wj,  vjf,  X%  =  [(X*)T,  ( * J")T,  Ct")T]T,  p  =  2  x  dim(.Vj) 


Initialization:  for  k  =  0,  set 

e+  =  E[0O] 

Xq  =  E[x0] 

r+0  =  E[(a-o-^)(x0-x+)t] 

Computation:  for  k  =  1,2,...  compute 
Parameter  estimate  time  update 
Parameter  covariance  time  update 


State  estimate  time  update 


State  covariance  time  update 
Output  estimate,  parameter  filter 


Output  estimate,  state  filter 


State  filter  gain  matrix 


Parameter  filter  gain  matrix 


State  estimate  measurement  update 
State  covariance  measurement  update 
Parameter  estimate  measurement  update 
Parameter  covariance  measurement  update 


=  E[(0o  -  §+)(0o  -  §+)T] 

Xq  +  =  E[.tg]  =  [(Xq  )t,  w,  D]t 

Ko  =  -  Sf+X*g  -  Xo'+)Tl  =  diag(rj0,  sv) 


6k  =  K-i 

JJ7  =  57+  +  E. 

e.k  O.k-l  r 


x“k-i  =  {*£}  •  *a-i  +  rfzftv  -  Y  } 

xkJ  =  /(*&.,.  xk-i,v  K-k~v 

^  =  TL  0«r^T 

^  =  Eto  -  Wa,t  -  ^)T 

wt  =  +  Y yjsfk,e~  -  Y 

©A,;  =  h(f(xk_ p  wa-1.  UFjfc_i,  W*,(,  *  -  1),  «a,  ®A.  W*,f,  A) 

4  =  EiLo  or  dk  =  ©A, o 

»  =  Ef=Xm%- 

w.a  =  Ef=0"Jc)%.  -  &)%;  -  %)T 
^y,k  =  Ef=0  “'c)(*A:r  -  Ww  -  »)T 


Ll  =  S, 


-1 

xy,k^y,k 


©a, A  =  Ejla«*  -  4)(©a,;  -  4)T 

A  =  Ef=o  -  §A  X®»,-  -  dk)T 
L°k  =  Vlk 

Xk  =  X~  +  Ljji(y*  -  %) 

Xtk  =  Slk  ~  LlSu(Lxk)T 
K  =  ®k  +  Lk(yk  -  dk) 
Kk=sik-Leksi^lf 


ically, 


The  joint/dual  SPKF  is  run  on  this  modified  model,  with  the 
“measured”  information  in  the  measurement  update  being 


yk  ~  OCV(zk)  -  Rik 
OCV(zk)  ~  yk  +  Rik 
Zk  =  OCV~l(yk  +  Rik)- 


By  measuring  the  cell  voltage  under  load,  the  cell  current,  and 
having  knowledge  of  R  (perhaps  through  0  from  the  dual/joint 
SPKF),  and  knowing  the  inverse  OC  V  function  for  the  cell  chem¬ 
istry,  one  can  compute  a  noisy  estimate  of  SOC,  Zk- 

To  combine  this  simple  model  with  the  dual/joint  SPKF, 
the  cell  model  being  used  (e.g.,  perhaps  the  ESC  model  in 
Section  7)  has  its  output  equation  augmented  with  SOC: 


h(xk,  Uk ,  vk,  Ok,  k)  = 


OC V(zk)  -  Rik  +  Gfk  +  Mhk 

Zk 


yk 

h 

While  the  “noise”  of  zk  (short-term  bias  due  to  hysteresis  ef¬ 
fects  and  polarization  filter  voltages  being  ignored)  prohibit  it 
from  being  used  as  the  primary  estimator  of  SOC,  its  expected 
long-term  behavior  in  a  dynamic  environment  is  accurate,  and 
maintains  the  accuracy  of  the  SOC  state  in  the  joint/dual  SPKF. 

7.  The  enhanced  self-correcting  cell  model 

In  order  to  examine  and  compare  performance  of  the  pro¬ 
posed  algorithms,  we  must  first  define  a  discrete-time  state- 
space  model  of  the  form  of  (1)  and  (2)  that  applies  to  battery 
cells.  Here,  we  briefly  review  the  “enhanced  self-correcting  cell 
model”  from  Refs.  [5,3].  This  model  includes  effects  due  to 
open-circuit- voltage,  internal  resistance,  voltage  time  constants, 
and  hysteresis.  For  the  purpose  of  example,  we  will  later  fit  pa- 
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SOC  as  a  function  of  time 


(a)  Time  (min.) 


Current  for  one  UDDS  cycle  (zoom) 


Fig.  2.  Plots  showing  SOC  vs.  time  and  rate  vs.  time  for  UDDS  cell  tests.  SOC  is  shown  in  (a)  and  rate  for  one  UDDS  cycle  is  shown  in  (b). 


rameter  values  to  this  model  structure  to  model  the  dynamics 
of  high-power  lithium-ion  polymer  battery  cells,  although  the 
structure  and  methods  presented  here  are  general. 

State-of-charge  is  captured  by  one  state  of  the  model.  This 
equation  is 

friiAT\  . 

Zk  =  Zk-l  -  I  — — —  I  IJfc-1, 

where  AT  represents  the  inter-sample  period  (in  seconds),  and 
C  represents  the  cell  capacity  (in  ampere-seconds). 

The  time  constants  of  the  cell  voltage  response  are  captured 
by  several  filter  states.  If  we  let  there  be  n  f  time  constants,  then 


fk  =  Affk-i  +  Bfik-\. 


The  matrix  Ar  e  ]R'!/X',/  may  be  a  diagonal  matrix  with  real¬ 
valued  entries.  If  so,  the  system  is  stable  if  all  entries  have  mag¬ 
nitude  less  than  one.  The  vector  B  f  e  R"Txl  may  simply  be  set 
to  n  f  “l”s.  The  value  of  n  f  and  the  entries  in  the  A  f  matrix  are 
chosen  as  part  of  the  system  identification  procedure  to  best  fit 
the  model  parameters  to  measured  cell  data. 

The  hysteresis  level  is  captured  by  a  single  state 


hk  —  exp 


rgik-iy  AT 
C 


hk-i 


+ 


1  —  exp 


rgik-W  AT 
C 


sgn(4-i), 


8.  Application  to  battery  management  systems 

8.1.  Cell  and  cell  test  description 

The  cells  used  in  this  paper  differ  electrochemically  from 
those  reported  in  previous  work.  We  refer  to  the  older  cells  as 
GEN3  cells,  and  to  the  newer  cells  as  G4  cells.  The  GEN3  cells 
are  high-power  (>20  C  capable)  7.5Ah  Mn  spinel/graphite  LiPB, 
and  the  G4  cells  are  very  high-power  (>30  C  capable)  5Ah  Mn 
spinel/blended-carbon  LiPB,  both  reported  in  Ref.  [22], 

In  order  to  compare  the  various  Kalman  filtering  methods’ 
abilities  to  estimate  SOC  and  SOH,  we  gathered  data  from  two 
prototype  LiPB  cells.  One  cell’s  data  was  used  to  tune  cell  model 
parameters,  and  the  second  cell’s  data  was  used  in  some  tests  to 
see  how  well  the  filters  generalized  to  slightly  different  dynamics 
than  expected.  For  the  tests,  we  used  a  Tenney  thermal  chamber 
set  at  25  °C  and  an  Arbin  BT2000  cell  cycler.  Each  channel  of 
the  Arbin  was  capable  of  20  A  current,  and  10  channels  were 
connected  in  parallel  to  achieve  currents  of  up  to  200  A.  The  cy¬ 
cler’s  voltage  measurement  accuracy  was  ±5  mV  and  its  current 
measurement  accuracy  was  ±200  mA. 

The  cell  test  we  use  here  comprised  a  sequence  of  1 8  (full) 
“urban  dynamometer  driving  schedule”  (UDDS)  cycles,  sepa¬ 
rated  by  15  A  discharge  pulses  and  5-min  rests,  and  spread  over 
the  90-10%  SOC  range.  The  SOC  as  a  function  of  time  is  plotted 
in  Fig.  2(a),  and  rate  as  a  function  of  time  for  one  of  the  UDDS 


where  y  is  the  hysteresis  rate  constant,  again  found  by  system 
identification. 

The  overall  model  state  is 

Xk  =  [fk,  hk,  H-]T- 

The  state  equation  for  the  model  is  formed  by  combining  all  of 
the  individual  equations,  above. 

The  output  equation  that  combines  the  state  values  to  predict 
cell  voltage  is 

yk  =  OCV(zfc)  ±  Gfk  -  Rik  +  Mhk, 

where  G  e  IlxV  is  a  vector  of  constants  that  blend  the  time- 
constant  states  together  in  the  output,  R  the  cell  resistance  (dif¬ 
ferent  values  may  be  used  for  dis/charge),  and  M  is  the  maximum 
hysteresis  level. 


OCV  as  a  function  of  SOC 


Fig.  3.  Plot  of  open-circuit-voltage  as  a  function  of  state-of-charge. 
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Voltage  for  one  UDDS  cvcle  (zoom)  True  SOC  and  voltage-based  estimate 


(a)  Time  (min.)  (b)  Time  (min.) 

Fig.  4.  Modeling  of  voltage  and  SOC:  (a)  ESC  modeling  of  cell  voltage  vs.  true  cell  voltage  for  one  UDDS  cycle  and  (b)  instantaneous  voltage-based  (non-SPKF-based) 
SOC  estimate  plotted  vs.  true  SOC. 


Estimation  error:  SPKF  ESC,  n=  2,  25°C  Estimation  error:  SPKF  ESC,  n=  2,  25°C 


0  60  120  180  240  300  360  420  480  540  600  0  60  120  180  240  300  360  420  480  540  600 

(a)  Time  (min.)  (c)  Time  (min.) 


Estimation  error:  SR-SPKF  ESC,  n=2,  25°C  Estimation  error:  SR-SPKF  ESC,  n  =  2,  25°C 


0  60  120  180  240  300  360  420  480  540  600  0  60  120  180  240  300  360  420  480  540  600 

(b)  Time  (min.)  (d)  Time  (min.) 


Fig.  5.  SOC  estimation  error  for  SPKF  vs.  SR-SPKF  with  correct  filter  initialization:  (a  and  b)  use  training  data  and  (c  and  d)  use  generalization  data. 


Table  8 

Comparison  of  SPKF  vs.  SR-SPKF  in  UDDS  test  results  predicting  SOC 

Correctly  initialized  Incorrectly  initialized 


RMS  error 

(%) 

Maximum 
error  (%) 

Bounds  error  (%) 

RMS  error  (%) 

Maximum  error  (%) 

Bounds  error  (%) 

Time  to  converge  (s) 

SPKF 

0.30 

1.51 

0.99 

1.44 

20.19 

2.74 

1138 

SR-SPKF 

0.30 

1.51 

0.98 

1.44 

20.19 

2.68 

1138 

SPKF  generalize 

0.95 

4.49 

0.35 

2.01 

23.15 

7.75 

2816 

SR-SPKF  generalize 

0.94 

4.49 

0.35 

2.01 

23.15 

7.76 

2816 
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Estimation  error:  SPKF  ESC,  n  =2,  25°C 


Estimation  error:  SPKF  ESC,  n=  2,  25°C 


(a)  Time  (min.)  (c)  Time  (min.) 


Estimation  error:  SR-SPKF  ESC,  n  =2,  25°C  Estimation  error:  SR-SPKF  ESC,  n  =2,  25°C 


(b)  Time  (min.)  (d)  Time  (min.) 


Fig.  6.  SOC  estimation  error  for  SPKF  vs.  SR-SPKF  with  incorrect  filter  initialization:  (a  and  b)  use  training  data  and  (c  and  d)  use  generalization  data. 


Table  9 

Comparison  of  joint  SPKF  vs.  dual  SPKF  vs.  dual  SR-SPKF  in  UDDS  test  results  predicting  SOC 


Correctly  initialized 

Incorrectly  initialized 

RMS  error  (%) 

Maximum 
error  (%) 

Bounds  error  (%) 

RMS  error  (%) 

Maximum  error  (%) 

Bounds  error  (%) 

Time  to 
converge  (s) 

Joint  SPKF 

0.29 

1.34 

1.31 

1.13 

19.83 

3.02 

1134 

Dual  SPKF 

0.32 

1.33 

0.21 

1.35 

19.83 

2.15 

1181 

Dual  SR-SPKF 

0.27 

1.33 

1.05 

1.26 

19.83 

2.74 

1138 

Joint  SPKF  generalize 

0.93 

4.40 

0.21 

2.11 

22.94 

8.02 

2839 

Dual  SPKF  generalize 

1.18 

5.43 

10.75 

2.02 

22.87 

12.30 

3950 

Dual  SR-SPKF  generalize 

0.90 

4.35 

0.23 

2.01 

22.87 

10.33 

3369 

cycles  is  plotted  in  Fig.  2(b).  We  see  that  SOC  increases  by  about 
5%  during  each  UDDS  cycle,  but  is  brought  down  slightly  less 
than  10%  during  each  discharge  between  cycles.  The  entire  op¬ 
erating  range  for  these  cells  ( 1 0-90%  SOC,  delineated  by  dashed 
lines  in  the  figure)  is  excited  during  the  cell  test. 

8.2.  Fitting  data  to  the  enhanced  self-correcting  model 

Data  collected  from  the  first  cell  was  used  to  identify  initial 
parameters  for  the  ESC  cell  model.  The  goal  is  to  have  the  cell 
model  output  resemble  the  cell  terminal  voltage  under  load  as 
closely  as  possible,  at  all  times,  when  the  cell  model  input  is 
equal  to  the  cell  current.  Model  fit  was  judged  by  comparing  root- 
mean-squared  (RMS)  estimation  error  (estimation  error  equals 
cell  voltage  minus  model  voltage)  over  the  portions  of  the  cell 
tests  where  SOC  was  between  5  and  95%.  Model  error  outside 


that  SOC  range  was  not  considered  as  the  HEV  pack  operation 
design  limits  are  10-90%  SOC.  Details  for  how  the  open-circuit- 
voltage  curve  was  generated  and  how  the  model  parameters  were 
fit  are  described  in  Ref.  [23].  In  particular,  the  open-circuit- 
voltage  as  a  function  of  state-of-charge  for  these  cells  is  plotted  in 
Fig.  3.  Values  were  fit  to  the  other  ESC-model  parameters,  with 
very  close  agreement  between  the  cell  model  voltage  prediction 
and  the  cell  true  voltage.  In  this  work,  the  model  employs  two 
low-pass  filter  states  (n  f  =  2),  a  nominal  capacity  of  5.0  Ah, 
and  an  inter-sample  interval  of  A  T  —  Is.4 


4  We  note  here  that  prior  work  with  the  third-generation  cells  used  n  /  =  4,  and 
that  this  work  using  fourth-generation  cells  uses  n  f  =  2.  There  are  competing 
objectives  here:  to  make  the  model  as  accurate  as  possible,  and  to  make  the  filter 
as  computationally  efficient  as  possible.  The  minimum  number  of  filter  states 
that  can  result  in  a  zero  dc-gain  is  two,  and  we  find  that  SPKF  is  enough  superior 
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Estimation  error:  Joint  SPKF  ESC,  nf=  2,  25  C 


Estimation  error:  Joint  SPKF  ESC,  n=  2,  25°C 


(a) 


(b) 


Time  (min.) 

Estimation  error:  Dual  SPKF  ESC,  nf  =2,  25°C 


Time  (min.) 

Estimation  error:  Dual  SR-SPKF  ESC,  n  =2,  25°C 


(d) 


(e) 


Time  (min.) 

Estimation  error:  Dual  SPKF  ESC,  n  =2,  25°C 


Time  (min.) 


Estimation  error:  Dual  SR-SPKF  ESC,  n  =2,  25°C 


(c) 


Time  (min.) 


(f) 


Time  (min.) 


Fig.  7.  SOC  estimation  error  for  joint  SPKF  vs.  dual  SPKF  vs.  dual  SR-SPKF  with  correct  initialization:  (a~c)  use  training  data  and  (d — f)  use  generalization  data. 


Some  modeling  results  are  shown  in  Fig.  4.  In  frame  (a),  the 
ESC  cell  voltage  prediction  is  compared  with  test  data.  Very 
close  agreement  is  observed  especially  during  the  dynamic  part 
of  the  test.  In  frame  (b),  an  instantaneous  voltage-based  SOC  es¬ 
timate,  calculated  as  z.k  —  OCV- 1  (yt  +  Rik)  (cf.  Section  6.1), 
is  compared  with  true  SOC.* * 5  This  comparison  is  made  to  show: 


to  EKF  that  it  gives  acceptable  estimates  with  a  poorer  model.  Therefore,  we 

have  chosen  to  use  nf  =  2  here. 

5  We  calculate  “true”  SOC  using  Coulomb  counting  from  the  Arbin  test  equip¬ 
ment.  Bias  and  noise  in  the  Arbin  current  sensor  will  cause  this  value  to  drift  from 
the  ideal  value,  but  over  the  relatively  short  duration  of  the  tests  and  given  the 


(1)  that  such  an  estimate  is  too  noisy  to  be  used  as  an  estimate 
of  SOC  by  itself  (and  therefore  we  need  to  use  more  advanced 
methods,  such  as  SPKF),  but  (2)  it  yields  average  behavior  in 
dynamic  tests  that  is  accurate,  and  so  is  useful  to  ensure  conver¬ 
gence  of  the  parameters  in  a  dual  or  joint  application.  Note  that 
there  is  no  point  in  low-pass  filtering  this  result,  as  the  delay  of 
such  a  filter  would  make  the  estimate  useless.  Rather,  the  tuning 


high  accuracy  of  the  Arbin  sensors  we  feel  that  this  is  a  reasonable  approximation 
to  true  SOC.  The  cost  of  using  such  high-accuracy  sensors  in  a  production  BMS 
is  prohibitive,  which  is  why  we  instead  investigate  using  intelligent  algorithms 
and  less-expensive  sensors. 
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Estimation  error:  Joint  SPKF  ESC,  r\f=  2,  25°C 


Estimation  error:  Joint  SPKF  ESC,  nf- 2,  25°C 


(d)  Time  (min.) 


Estimation  error:  Dual  SPKF  ESC,  nf=  2,  25°C 


Time  (min.) 


Estimation  error:  Dual  SPKF  ESC,  n=  2,  25°C 


(e)  Time  (min.) 


Estimation  error:  Dual  SR-SPKF  ESC,  nf=  2,  25°C 


Estimation  error:  Dual  SR-SPKF  ESC,  nf= 2,  25°C 


Fig.  8.  SOC  estimation  error  for  joint  SPKF  vs.  dual  SPKF  vs.  dual  SR-SPKF  with  incorrect  initialization:  (a-c)  use  training  data  and  (d-f)  use  generalization  data. 


parameters  of  the  SPKF  are  used  to  adjust  the  “belief”  that  the 
filter  has  in  the  accuracy  of  z.k  when  adapting  its  own  internal 
state. 

8.3.  Examples  of  SPKF  and  SR-SPKF 

In  the  companion  to  this  paper,  we  showed  SPKF  SOC  esti¬ 
mation  results  compared  against  EKF  results  for  th  GEN3  cell. 
The  conclusion  was  that  the  SPKF  was  better  in  all  cases,  so  we 
do  not  include  EKF  results  here.  Rather,  we  begin  by  present¬ 
ing  SPKF  results  for  the  G4  cell.  For  brevity,  we  include  only 
room-temperature  results  here. 


Fig.  5(a)  shows  SOC  estimation  error  and  predicted  error 
bounds  for  the  training  cell,  where  the  SPKF  state  was  correctly 
initialized.6  In  particular,  the  SOC  state  of  the  SPKF  was  ini¬ 
tialized  to  the  correct  value  of  100%.  Frame  (b)  shows  the  same 
result  using  SR-SPKF.  The  results  are  nearly  identical,  as  ex¬ 
pected.  RMS  SOC  estimation  error  was  0.3%  and  maximum 


6  The  error  bounds  are  computed  as  plus/minus  three  times  the  square-root  of 
the  covariance  matrix  diagonal  element  corresponding  to  the  SOC  state.  These 
are  referred  to  as  "three-sigma”  bounds,  and  if  all  densities  are  Gaussian,  the 
bounds  should  correctly  encompass  the  true  value  of  SOC  99.7%  of  the  time. 
Section  8.5  discusses  the  validity  of  this  assumption  in  more  detail. 
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Instantaneous  resistance  and  capacity  estimates 


Time  (min.) 

Fig.  9.  Instantaneous  resistance  and  capacity  estimates  by  joint  SPKF  (dual  filter 
results  similar). 

SOC  error  was  1.51%.  The  estimation  error  bounds  correctly 
included  the  true  SOC  value  for  all  but  about  1%  of  the  test  run. 
The  bounds  were  incorrect  only  at  the  very  end  of  the  test  where 
true  SOC  was  below  10%  and  out  of  the  expected  operating 
range  of  the  cell — thus  the  bounds  correctly  included  the  true 
SOC  for  all  SOCs  of  importance. 

Frames  (c)  and  (d)  show  the  same  results  when  the  SPKF 
was  applied  to  data  from  the  second  cell — the  one  not  used  to 


SOC  error:  Training  data  and  good  initialization 


(a)  Standard  deviations  away  from  zero 


SOC  error:  Testing  data  and  good  initialization 


(c)  Standard  deviations  away  from  zero 


train  the  model  parameters.  We  see  slightly  poorer  performance, 
as  would  be  expected.  Note  that  post-test  analysis  indicates  that 
this  cell  may  have  been  overcharged  to  about  104%  before  the 
test  began  -  based  on  OC  V  data  -  which  explains  why  the  SPKF 
gave  an  initial  4%  error  when  initialized  to  100%.  The  filter  was 
still  able  to  converge  to  the  true  SOC. 

Fig.  6  shows  results  parallel  to  those  of  Fig.  5  for  the  case 
where  the  filter  SOC  state  was  intentionally  initialized  with  an 
incorrect  value  (80%  rather  than  100%)  to  demonstrate  the  con¬ 
vergence  properties  of  SPKF.  We  see  relatively  fast  convergence 
to  within  5%  SOC  error,  about  20-min  convergence  to  2%  error 
for  the  training  cell  and  about  45-min  convergence  to  2%  error 
for  the  generalization  cell.  Convergence  time  may  be  adjusted 
by  varying  Ew  and  Uv,  and  could  have  been  made  faster  than 
the  results  shown.  The  tradeoff  would  have  been  poorer  SOC 
bounds  estimate. 

Table  8  summarizes  the  results  of  these  tests.  The  “bounds 
error”  column  shows  the  percentage  of  time  the  SOC  estimation 
bounds  did  not  encompass  the  true  SOC.  Note  that  this  was 
always  at  very  low  and  high  SOC  -  out  of  the  operating  range 
of  the  cell  -  so  is  not  of  concern  here.  The  “time  to  converge” 
column  for  the  generalization  tests  show  how  many  iterations 
(seconds)  were  required  for  the  SOC  error  to  converge  to  less 
than  2%. 


SOC  error:  Training  data  and  bad  initialization 


SOC  error:  Testing  data  and  bad  initialization 


(d)  Standard  deviations  away  from  zero 


Fig.  10.  Empirical  probability-density  functions  of  SOC  estimation  error  normalized  by  standard  deviation:  (a)  training  data  and  correct  initialization,  (b)  training  data 
and  incorrect  initialization,  (c)  generalization  data  and  correct  initialization,  and  (d)  generalization  data  and  incorrect  initialization.  The  standard-normal  distribution 
is  overlaid  in  each  case. 
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8.4.  Examples  of  joint  and  dual  ID 

We  conducted  tests  to  determine  the  SOC  estimation  perfor¬ 
mance  of  the  joint/dual  SPKF  filters  as  well.  Figs.  7  and  8  show 
results  parallel  to  those  of  Figs.  5  and  6. 

In  Fig.  7,  results  are  shown  when  the  filter  state  was  correctly 
initialized;  in  Fig.  8,  the  SOC  state  was  intentionally  initialized  to 
an  incorrect  value  (80%  rather  than  100%)  to  show  convergence 
properties.  In  frames  (a)-(c),  results  using  the  cell  from  which 
parameters  were  originally  trained  are  shown;  in  frames  (d)- 
(f),  generalization  results  are  shown.  Note  that  the  joint/dual 
filters  adapt  their  parameters,  so  that  over  time  they  learn  the 
dynamics  of  any  cell  they  are  modeling,  and  SOC  estimation 
error  improves. 

Table  9  summarizes  results  for  these  tests.  We  note  that  they 
are  very  similar  to  the  SOC  results  for  the  plain  SPKF  filter.  Over 
a  (much)  longer  test  we  would  expect  some  improvement  in  the 
joint/dual  results  as  the  filters  adapt  to  changing  cell  dynamics. 

In  particular,  Fig.  9  shows  the  1-Hz  resistance  and  capacity 
parameters  adapting  over  the  course  of  the  test.  The  change  in  re¬ 
sistance  is  valid — these  cells  have  much  higher  resistance  at  low 
SOC  than  at  moderate  SOC.  The  variation  in  the  capacity  esti¬ 
mate  between  5.02  and  5.04  Ah  is  a  little  less  certain.  It  is  likely 
that  the  capacity  state  is  exhibiting  some  short-term  variation 
to  mask  some  other  modeling  error.  We  have  observed  that  the 
long-term  behavior  of  the  capacity  state  is  stable,  however,  and 
is  able  to  track  changes  in  the  true  cell  capacity.  The  adaptation 
may  be  slowed  down  to  minimize  the  masking  effect  shown  here 
by  lowering  the  fictitious  noise  covariance  in  Er  corresponding 
to  the  capacity  state. 

8.5.  Analysis  of  SOC  estimation  error 

Before  concluding  this  paper,  we  present  some  results  ana¬ 
lyzing  the  SOC  estimation  error.  All  of  the  Kalman  filter  vari¬ 
ants  make  the  assumption  that  the  noises  affecting  the  system 
of  concern  and  the  measurements  made  are  Gaussian  (normal) 
random  variables,  and  that  the  state  estimation  errors  are  also 
Gaussian.  In  Fig.  10,  we  plot  summary  histograms  of  SOC  es¬ 
timation  error  accumulated  through  the  various  test  runs  using 
SPKF,  SR-SPKF,  joint  SPKF,  dual  SPKF,  and  dual  SR-SPKF.  In 
every  case,  each  instantaneous  SOC  estimation  error  was  divided 
by  its  corresponding  one-sigma  error  bound  before  computing 
the  histogram,  and  then  the  histogram  was  normalized  to  have 
unit  area  to  form  an  empirical  probability  density  function  (PDF) 
of  the  SOC  estimation  error.  If  all  KF  assumptions  were  being 
met,  this  distribution  should  match  a  unit-variance  standard  nor¬ 
mal  PDF.  We  see  that  this  assumption  is  fairly  close  to  being  met 
when  using  the  training  data,  but  less  well  met  when  using  the 
generalization  data.  In  particular,  we  note  that  three-sigma  error 
bounds  (as  reported  elsewhere  in  this  paper)  are  theoretically 
accurate  99.7%  of  the  time,  but  we  see  from  the  histograms  that 
four-sigma  bounds  would  give  a  little  more  safety  in  an  imple¬ 
mentation.  Also,  since  we  now  see  that  the  KF  assumptions  are 
not  being  met,  doubt  is  cast  on  how  much  better  we  might  be 
able  to  do  using  KF  techniques.  Most  likely,  one  would  need  to 
use  a  particle  filter  -  which  does  not  assume  Gaussian  RVs  -  to 


do  better,  but  the  added  computational  complexity  might  not  be 
tolerable  in  a  commercial  application. 

9.  Conclusions 

This  paper  concludes  a  two-part  series  discussing  the  ap¬ 
plication  of  sigma-point  Kalman  filters  to  battery  management 
algorithms.  In  the  first  paper,  we  introduced  the  general  prob¬ 
abilistic  inference  solution  to  optimal  estimation,  and  derived 
the  KF,  EKF,  and  SPKF  from  this  solution  using  different  sets 
of  assumptions.  The  SPKF  was  shown  to  be  theoretically  more 
precise  than  EKF;  testing  with  real  cell  data  supported  this  anal¬ 
ysis. 

This  paper  showed  how  SPKF  could  be  very  closely  approx¬ 
imated  by  SR-SPKF,  which  gives  speed  advantages  in  the  case 
of  a  linear  state  equation,  such  as  when  estimating  parameters 
rather  than  states.  Simultaneous  state  and  parameter  estimation 
was  also  introduced  via  joint  and  dual  SPKF  methods.  Results 
were  presented  for  a  prototype  LiPB  HEV  cell  that  demonstrate 
that  these  methods  work  very  well.  For  example,  typical  SOC 
estimation  errors  of  less  than  1%  are  reported,  both  using  train¬ 
ing  and  testing  data,  with  near  perfect  error  bounds — these  were 
in  fact  perfect  over  the  10-90%  expected  SOC  operating  range 
of  the  cell.  Parameters  such  as  cell  resistance  and  SOC  can  also 
be  simultaneously  estimated. 

The  various  KF  methods  are  derived  assuming  that  the  prob¬ 
ability  density  functions  of  sensor  and  process  noises  are  Gaus¬ 
sian.  We  have  seen  by  example  that  this  assumption  is  not  strictly 
adhered  to  in  this  application  (nor  can  it  be  exactly).  This  limits 
the  ability  of  the  KF  methods  to  estimate  states  and  parameters, 
and  indicates  that  different  approaches  might  need  to  be  taken  if 
even  greater  accuracy  is  needed.  However,  the  SPKF  does  very 
well  even  though  this  assumption  is  not  met  exactly. 

For  further  reading,  we  have  shown  elsewhere  that  the  state 
and  parameters  may  be  used  to  very  precisely  estimate  dynamic 
available  power,  and  to  compute  which  cells  must  be  equalized 
[6,24]. 
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